{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Combustion Model and Ignition Delay Demo\n",
    "\n",
    "Written by A. Mark Payne and Alon Grinberg Dana for presentation at ICCK 2017"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## User input:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Fuel = 'CCO'\n",
    "Equivalence_Ratio = 1.0\n",
    "\n",
    "Temperature = 1500.0  # (K)\n",
    "Pressure = 1.0        # (atm)\n",
    "simtime = 2           # (ms)\n",
    "topSA = 10            # number of top sensitive reactions and thermo to display\n",
    "\n",
    "rmgpy_path = '../rmg.py' # Change to your rmg.py path\n",
    "\n",
    "\n",
    "from IPython.display import display, Image\n",
    "from rmgpy.molecule import Molecule\n",
    "Fuel_Molecule = Molecule().fromSMILES(Fuel)\n",
    "print(\"The fuel molecule is:\");display(Fuel_Molecule)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RMG's input file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Fuel_Molecule = Molecule().fromSMILES(Fuel)\n",
    "nC = int(Fuel_Molecule.getNumAtoms('C'))\n",
    "nH = int(Fuel_Molecule.getNumAtoms('H'))\n",
    "nO = int(Fuel_Molecule.getNumAtoms('O'))\n",
    "\n",
    "A = str(Equivalence_Ratio/(nC+(nH/4.0)-(nO/2.0)))\n",
    "\n",
    "Input_file = '''\n",
    "# Data sources\n",
    "database(\n",
    "    thermoLibraries = ['BurkeH2O2','primaryThermoLibrary','thermo_DFT_CCSDTF12_BAC','DFT_QCI_thermo','FFCM1(-)','JetSurF2.0'],\n",
    "    reactionLibraries = ['BurkeH2O2inN2','FFCM1(-)','JetSurF2.0'],\n",
    "    seedMechanisms = [],\n",
    "    kineticsDepositories = ['training'],\n",
    "    kineticsFamilies = 'default',\n",
    "    kineticsEstimator = 'rate rules',\n",
    ")\n",
    "\n",
    "# List of species\n",
    "species(\n",
    "    label='fuel',\n",
    "    reactive=True,\n",
    "    structure=SMILES('''+\"'\"+Fuel+\"'\"+'''),\n",
    ")\n",
    "\n",
    "species(\n",
    "    label='O2',\n",
    "    reactive=True,\n",
    "    structure=SMILES('[O][O]'),\n",
    ")\n",
    "\n",
    "species(\n",
    "    label='N2',\n",
    "    reactive=False,\n",
    "    structure=SMILES('N#N'),\n",
    ")\n",
    "\n",
    "species(\n",
    "    label='OH',\n",
    "    reactive=True,\n",
    "    structure=SMILES('[OH]'),\n",
    ")\n",
    "\n",
    "# Reaction system\n",
    "simpleReactor(\n",
    "    temperature=('''+str(Temperature)+''','K'),\n",
    "    pressure=('''+str(Pressure)+''','atm'),\n",
    "    initialMoleFractions={\n",
    "        'fuel': '''+A+''',\n",
    "        'O2': 1,\n",
    "        'N2': 3.76,\n",
    "    },\n",
    "    terminationTime=(0.001,'s'),\n",
    "    sensitivity=['OH'],\n",
    "    sensitivityThreshold=0.01,\n",
    ")\n",
    "\n",
    "simulator(\n",
    "    atol=1e-16,\n",
    "    rtol=1e-8,\n",
    "    sens_atol=1e-6,\n",
    "    sens_rtol=1e-4,\n",
    ")\n",
    "\n",
    "model(\n",
    "    toleranceKeepInEdge=0,\n",
    "    toleranceMoveToCore=0.05,\n",
    "    toleranceInterruptSimulation=0.05,\n",
    "    maximumEdgeSpecies=300000\n",
    ")\n",
    "\n",
    "#pressureDependence(\n",
    "#        method='modified strong collision',\n",
    "#        maximumGrainSize=(0.5,'kcal/mol'),\n",
    "#        minimumNumberOfGrains=250,\n",
    "#        temperatures=(298,2500,'K',10),\n",
    "#        pressures=(0.5,3,'bar',5),\n",
    "#        interpolation=('Chebyshev', 6, 4),\n",
    "#        maximumAtoms=16,\n",
    "#)\n",
    "\n",
    "options(\n",
    "    units='si',\n",
    "    generateOutputHTML=True,\n",
    "    generatePlots=False,\n",
    "    saveEdgeSpecies=False,\n",
    "    saveSimulationProfiles=True,\n",
    "    saveRestartPeriod=None,\n",
    ")\n",
    "\n",
    "generatedSpeciesConstraints(\n",
    "    allowed=['input species','seed mechanisms','reaction libraries'],\n",
    "    maximumCarbonAtoms=5,\n",
    "    maximumOxygenAtoms=2,\n",
    "    maximumNitrogenAtoms=0,\n",
    "    maximumSiliconAtoms=0,\n",
    "    maximumSulfurAtoms=0,\n",
    "    maximumHeavyAtoms=6,\n",
    "    maximumRadicalElectrons=2,\n",
    "    allowSingletO2=False,\n",
    ")\n",
    "\n",
    "\n",
    "'''\n",
    "\n",
    "import os\n",
    "if not os.path.exists('RMG'):\n",
    "    os.mkdir('RMG')\n",
    "os.system('rm -r RMG/*')\n",
    "with open('RMG/Demo.py','w') as RMG_Input_File:\n",
    "    RMG_Input_File.write(Input_file)\n",
    "print(\"Created RMG input file\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run RMG:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "os.system('python {0} RMG/Demo.py'.format(rmgpy_path))         \n",
    "\n",
    "print(\"RMG Simulation Completed. Summary of log file:\\n\")\n",
    "RMG_log = open('RMG/RMG.log','r').readlines()\n",
    "lines = [x for x in RMG_log[-13:-1] if x != \"\\n\"]\n",
    "for line in lines: print(line)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run the generated model (using RMG's Cantera functions):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from rmgpy.chemkin import *\n",
    "from rmgpy.tools.canteraModel import *\n",
    "from rmgpy.species import Species\n",
    "import time\n",
    "\n",
    "path = \"RMG/chemkin/\"\n",
    "speciesList, reactionList = loadChemkinFile(path+'chem_annotated.inp',\n",
    "                                            path+'species_dictionary.txt',\n",
    "                                            path+'tran.dat')\n",
    "\n",
    "nC = int(Fuel_Molecule.getNumAtoms('C'))\n",
    "nH = int(Fuel_Molecule.getNumAtoms('H'))\n",
    "nO = int(Fuel_Molecule.getNumAtoms('O'))\n",
    "phi = Equivalence_Ratio\n",
    "FuelStoich = phi/(nC+(nH/4.0)-(nO/2.0))\n",
    "\n",
    "Fuel_Species=Species().fromSMILES(Fuel)\n",
    "O2_Species=Species().fromSMILES('[O][O]')\n",
    "N2_Species=Species().fromSMILES('N#N')\n",
    "OH_Species=Species().fromSMILES('[OH]')\n",
    "species_dict = getRMGSpeciesFromUserSpecies([Fuel_Species,O2_Species,N2_Species,OH_Species], speciesList)\n",
    "\n",
    "reactorTypeList = ['IdealGasReactor']\n",
    "reactionTimeList = ([simtime], 'ms')\n",
    "\n",
    "molFracList=[{species_dict[Fuel_Species]: FuelStoich,\n",
    "              species_dict[O2_Species]: 1,\n",
    "              species_dict[N2_Species]: 3.76}]\n",
    "Tlist = ([Temperature],'K')\n",
    "Plist = ([Pressure],'atm')\n",
    "\n",
    "job = Cantera(speciesList=speciesList, reactionList=reactionList, outputDirectory='')\n",
    "job.loadChemkinModel(path+'chem_annotated.inp',transportFile=path+'tran.dat')\n",
    "job.generateConditions(reactorTypeList, reactionTimeList, molFracList, Tlist, Plist)\n",
    "\n",
    "alldata = job.simulate()\n",
    "print(\"Done.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "############### Settings ###############\n",
    "fsize = (8,4) # Change to make the figure fit on your screen\n",
    "########################################\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from rmgpy.tools import plot as rmg_plot\n",
    "from operator import itemgetter\n",
    "from rmgpy.tools.sensitivity import runSensitivity\n",
    "%matplotlib notebook\n",
    "\n",
    "times = alldata[0][0].data\n",
    "pressures = alldata[0][1][1].data\n",
    "\n",
    "dpdt = (pressures[1:] - pressures[:-1]) / (times[1:] - times[:-1])\n",
    "idi = next(i for i,d in enumerate(dpdt) if d==max(dpdt))\n",
    "ign_delay_time = times[idi]\n",
    "\n",
    "for spc in xrange(len(alldata[0][1][:])):\n",
    "    if alldata[0][1][spc].label == str(species_dict[Fuel_Species]):\n",
    "        Fuel_idx = spc\n",
    "    if alldata[0][1][spc].label == str(species_dict[OH_Species]):\n",
    "        OH_idx = spc\n",
    "\n",
    "for i in range(len(alldata[0][1][Fuel_idx].data)):\n",
    "    if alldata[0][1][Fuel_idx].data[i]<0.001:\n",
    "        Fuel_Depletion_Time = times[i]\n",
    "        break\n",
    "\n",
    "files = os.listdir('RMG/solver')\n",
    "\n",
    "sensitivity_file = str(filter(lambda x: ('sensitivity' in x) and ('.csv' in x),files)[0])\n",
    "SA_time, SA_data = rmg_plot.parseCSVData('RMG/solver/'+sensitivity_file)\n",
    "\n",
    "time_error = 1\n",
    "\n",
    "for i in range(len(SA_time.data)):\n",
    "    if abs(SA_time.data[i]-ign_delay_time)<time_error:\n",
    "        ign_delay_idx = i\n",
    "        time_error = abs(SA_time.data[i]-ign_delay_time)\n",
    "\n",
    "Gidx = 0\n",
    "for i in xrange(len(SA_data[:])):\n",
    "    if \"G\" in SA_data[i].label:\n",
    "        if not Gidx:\n",
    "            Gidx = i\n",
    "        SA_data[i].label = SA_data[i].label.split('G')[1][1:-1]\n",
    "    else:\n",
    "        SA_data[i].label = SA_data[i].label.split(' ')[1]\n",
    "\n",
    "rank1 = []\n",
    "for n in xrange(Gidx):\n",
    "    rank1.append(abs(SA_data[n].data[ign_delay_idx]))               # list of max SA range for each rxn\n",
    "num1 = np.linspace(0,len(rank1)-1,len(rank1))        # will be used to get the order of reactions by rank1\n",
    "num1 = zip(rank1,num1)\n",
    "num1 = sorted(num1, key=itemgetter(0),reverse=True)\n",
    "SA_k_data = []\n",
    "SA_k_label = []\n",
    "for i in xrange(min(topSA, Gidx)):\n",
    "    SA_k_data.append(SA_data[int(num1[i][1])].data[ign_delay_idx])     # make sorted lists size topSA of SA values and rxns labels\n",
    "    SA_k_label.append(SA_data[int(num1[i][1])].label)\n",
    "        \n",
    "rank2 = []\n",
    "for n in xrange(len(SA_data)-Gidx):\n",
    "    rank2.append(abs(SA_data[n+Gidx].data[ign_delay_idx]))     # list of max SA range for each rxn\n",
    "num2 = np.linspace(0,len(rank2)-1,len(rank2))        # will be used to get the order of reactions by rank1\n",
    "num2 = zip(rank2,num2)\n",
    "num2 = sorted(num2, key=itemgetter(0),reverse=True)\n",
    "SA_G_data = []\n",
    "SA_G_label = []\n",
    "for i in xrange(min(topSA, len(SA_data)-Gidx)):\n",
    "    SA_G_data.append(SA_data[int(num2[i][1])+Gidx].data[ign_delay_idx])     # make sorted lists size topSA of SA values and rxns labels\n",
    "    SA_G_label.append(SA_data[int(num2[i][1])+Gidx].label)\n",
    "        \n",
    "plt.rcParams['axes.labelsize'] = 18\n",
    "plt.rcParams['xtick.labelsize'] = 12\n",
    "plt.rcParams['ytick.labelsize'] = 12\n",
    "plt.rcParams['figure.autolayout'] = True\n",
    "\n",
    "plt.style.use('ggplot')\n",
    "plt.style.use('seaborn-pastel')\n",
    "\n",
    "fig = plt.figure(figsize=fsize)\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "plt.plot(alldata[0][0].data*1000, alldata[0][1][Fuel_idx].data,'-o')\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('$Y_{Fuel}$')\n",
    "plt.title('Fuel profile')\n",
    "plt.xlim([0,2000*ign_delay_time])\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.plot(alldata[0][0].data*1000, alldata[0][1][OH_idx].data,'-o')\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('$Y_{OH}$')\n",
    "plt.title('OH profile')\n",
    "plt.xlim([0,2000*ign_delay_time])\n",
    "plt.arrow(0, alldata[0][1][OH_idx].data[idi], ign_delay_time*1000, 0, width=0.0001, head_width=0.0005, head_length=0.001, length_includes_head=True, color='r', shape='full')\n",
    "plt.annotate(r'$Ignition Delay: \\tau_{ign}$', xy=(0,0), xytext=(0, alldata[0][1][OH_idx].data[idi]+0.0005), fontsize=10);\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=fsize)\n",
    "plt.barh(np.arange(min(topSA, Gidx)), SA_k_data, 1/1.5, color=\"blue\")\n",
    "plt.gca().invert_yaxis()\n",
    "plt.xlabel(r'Sensitivity: $\\frac{\\partial\\:\\ln{[OH]}}{\\partial\\:\\ln{k}}$');\n",
    "plt.rcParams.update({'axes.labelsize': 20})\n",
    "plt.yticks(np.arange(min(topSA, Gidx)),SA_k_label)\n",
    "plt.title(\"[OH] sensitivity to kinetics\")\n",
    "\n",
    "fig = plt.figure(figsize=fsize)\n",
    "plt.barh(np.arange(min(topSA, len(SA_data)-Gidx)), SA_G_data, 1/1.5, color=\"blue\")\n",
    "plt.gca().invert_yaxis()\n",
    "plt.xlabel(r'Sensitivity: $\\frac{\\partial\\:\\ln{[OH]}}{\\partial\\:G_i}$ $[mol/kcal]$');\n",
    "plt.rcParams.update({'axes.labelsize': 20})\n",
    "plt.yticks(np.arange(min(topSA, len(SA_data)-Gidx)),SA_G_label)\n",
    "plt.title(\"[OH] sensitivity to thermo\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gas = ct.Solution('RMG/cantera/chem.cti')\n",
    "comp = str(species_dict[Fuel_Species])+\":\"+str(FuelStoich)+\",\"+str(species_dict[O2_Species])+\":1,\"+str(species_dict[N2_Species])+\":3.76\"\n",
    "gas.TPX = Temperature, Pressure, comp\n",
    "reactor = ct.IdealGasConstPressureReactor(gas)\n",
    "network = ct.ReactorNet([reactor])\n",
    "network.advance(ign_delay_time)\n",
    "ROP_C = ct.ReactionPathDiagram(gas, 'C')\n",
    "\n",
    "from PIL import Image as PILimg\n",
    "ROP1 = plt.subplot(1,1,1)\n",
    "dot_file = 'RMG/cantera/rxnpathC.dot'\n",
    "img_file = 'RMG/cantera/rxnpathC.png'\n",
    "ROP_C.title = 'Reaction path diagram following C'\n",
    "ROP_C.threshold = 0.01\n",
    "ROP_C.label_threshold = 0.01\n",
    "ROP_C.show_details = True\n",
    "ROP_C.write_dot(dot_file)                                              # write dot file\n",
    "os.system('dot {0} -Tpng -o{1} -Gdpi=300'.format(dot_file, img_file))   # write png file\n",
    "fullpath = os.getcwd() + '/' + img_file\n",
    "Image(fullpath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
